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Vacancy drag: A model for 
kinetics of the interface on the 
nanoscale 

4.1 Introduction and summary 

In the previous two chapters, I have investigated interfacial motion in the normal grain 
growth and in the presence of secondary-phase particles from a macroscopic perspective. 
In this chapter dynamics of interface with more concerns about its atomistic structure 
are considered. 

The geometric properties of interfaces in polycrystalline materials are usually im¬ 
posed by production history of the grains, e.g. solidification, cold work, heat treatment, 
etc. Thereafter, the grains have very little chance to move and adjust to the constrained 
surroundings. Hence the atoms located at the boundaries barely fit to the space be¬ 
tween the grains and there are localized free volumes at the interfaces. High diffusivity, 
high electrical resistivity, high chemical activity and toughness are general properties 
related to the irregularity and distortion at the grain boundaries. 

Corresponding to the free volumes at the interface, there are local elastic strains 
which may restrict the motion of the interface. One mechanism to relax the mechanical 
distortion is adsorbing (or emitting) mobile defects such as vacancies. A vacancy is a 
structural point defect of removing one atom from the lattice which carries small amount 
of strain. It is also the key element of self-diffusion. 

Solute atoms are also able to interact with interfaces and influence their migration 1 
either chemically or mechanically. Diffusion-induced interface motion is reported for 
some metallic polycrystalline materials in which accumulation of solute atoms at the 
interface motivates its migration. For the first time, Hillert and Purdy [82] reported this 
effect for the Fe-Zn system. The mechanism of this effect is widely studied. Hillert [83] 
and Penrose [84] discussed an elastic mechanism related to the lattice misfit between 
solute and solvent in which, the mechanical driving force is proportional to squared of 
concentration difference between the two sides of the interface, while Balluffi and Cahn 

1 Vacancies can be imagined as a secondary element (solute) within a pure substance. 
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[85] proposed that the difference in diffusion coefficients at the interface is responsible 
for diffusion-induced migration of the interface. 

Another study by Schmitz et al. [86] has shown that in a diffusion-couple, the order 
of layering may cause stresses which influence the diffusion rates: they have deposited 
nanosized Al/Cu/Al and Cu/Al/Cu layers on a finger tip of about 50 nm in diame¬ 
ter (curved geometry) and observed a significant imbalance of vacancy concentration 
existing in different layers. Owing to the geometry and small size of the system the 
internal stresses have been produced which lead to an asymmetric diffusion process in 
the system [86, 87]. 

In this chapter a pure nanocrystalline substance in which vacancies will be the only 
present point defects is considered. In Ns materials, the population of interfaces is very 
high and their free volumes no longer can be neglected. The interdependency between 
diffusion of vacancies and the stress build around the interfaces is taken into account 
and a phase-field model for the phenomenon of nanograin growth is constructed. The 
description of this model was first proposed by Steinbach et al. [88]. The effect of 
internal stress on the diffusion in solids, indeed, has been developed by Larche and 
Calm [34] and Khachaturian [35] and later has been used in different studies such as 
nucleation and growth of precipitates [89, 90]. 

The mechanical stresses can modify the diffusion mechanism in two ways: i) if 
the species create lattice distortions or ii) if the elastic constants are concentration- 
dependent, the diffusion can be driven by stress/strain gradients such that the total 
elastic strain is reduced. This is a delicate interdependency between concentration and 
stress fields. In the first case, the “mobile eigenstrains” compensate each other or any 
other source of distortion by means of diffusion. In a dynamic system consisting of 
many sinks and sources of distortion, this can lead to a complex trajectory of atoms 
and microstructural changes. In the second case, materials respond to the mechanical 
stress by local softening and hardening which is controlled by diffusion of suited species. 

In the following parts, I first give an overview about properties of Ns polycrystalline 
materials. Before introducing the methods, I briefly summarize the thermodynamic and 
mechanical characteristics of vacancies which seems necessary for our future reference. 
The results of the model are presented and compared with the background of other 
relevant models and experimental studies. In the present model, the dynamic of the 
interface is tuned by natural cooperation of the curvature, mechanical energy due to 
free volumes, and diffusion of vacancies. The curvature-driven motion of the interfaces 
is hindered by mechanical distortion. In order to unlock the interfaces, the interfa¬ 
cial distortion should be relaxed by diffusion of vacancies. Thereupon, there are two 
competitive kinetic processes, namely, the curvature-driven motion and the diffusion of 
vacancies which determine the kinetics of growth. Diffusion-controlled growth is found 
to be dominant on the nanoscale, specially at low temperatures. 


4.2 Nanocrystalline materials 

The object of study in this chapter is the contribution of diffusion (vacancies) to the 
grain growth in the presence of internal stress gradients. This is directly related to 
the nanocrystallined bulk materials which usually contain a large amount of defects 
and distortion. In this section, I briefly review the definition, properties and relevant 
findings about these special materials. 
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In the recent decades, Ns materials with an average grain size of less than 0.1 /mi 
have been considered as new possibilities for achieving better mechanical and functional 
properties [91, 92, 93, 94, 95]. The interfacial volume in Ns materials can reach up to 
50 % of the total volume and the density at the interfaces may drop to 70 %. A typical 
5nm-radius metallic grain contains only about 10,000 atoms. This grain may show a 
very different behaviour from a conventional grain with 1 /mi radius. For instance, 
the bulk of the nanograins is so small that there would be no source of dislocations to 
activate the slip systems under stress. Instead of that, the grain may compensate the 
external stress via rotation to a certain extent. The fabrication of Ns materials causes 
some side effects on the behaviour of grain boundaries as well [95, 96]. For example, 
mechanical procedures, such as mechanical alloying, induce a considerable amount of 
structural defects in atomistic and micro- levels which can influence following processes 
such as sintering or shaping. 

The large amount of point defects and interfaces provides short-circuits which en¬ 
hance the overall self-diffusion process. At room temperature, the self-diffusion for Ns 
copper is reported to be about 10 10 and 10 4 times larger than bulk and grain bound¬ 
ary diffusion, respectively. Increasing the temperature, the differences will be slightly 
less, but still significant [94], It is likely that the large number of triple junctions also 
play an important role in the self-diffusion mechanism. Solutes, however, may diffuse 
differently. Since there is a large amount of interfaces, the adsorption capacity is sig¬ 
nificantly high (high solute solubility) and the solute atoms may be easily trapped. 
However, once the adsorption sites are occupied at high concentration of solutes, the 
solute diffusivity may increase. It is even observed that once the concentration at the 
interfaces are high, the diffusion of solutes may propel the interface [82], This is called 
diffusion-induced motion. 

Grain size in polycrystalline materials strongly influences their mechanical prop¬ 
erties. The smaller the grain size, the stronger the materials is expected. At the 
nanoscale, however, this is accounted to be inverse below a certain size. This effect, 
known as inverse Hall-Petch effect, is attributed to sliding at the grain boundary in¬ 
stead of slip systems and also the lack of dislocations. Disordered arrangement of atoms 
at the interface may also change the heat transfer and heat capacity of these materials. 
It is reported that the heat capacity of some nanocrystallined materials, such as pal¬ 
ladium and copper, is more than corresponding conventional bulk material. For more 
information the reader is referred to [94, 95]. 

To apply nanonraterials, the stability of microstructure at desired temperatures 
must be warranted. The grain growth in nanocrystallines, however, differs from cor¬ 
responding coarse-grained bodies [6, 97, 98, 99]. The obscure and unpredictable mi- 
crostructural evolution in Ns materials requires solid studies to uncover the origins of 
nanograin growth. 

Different studies report a broad spectrum of observations concerning nanograin 
growth. Isothermal annealing experiments revel a slow linear growth regime at low 
temperatures that changes to normal grain growth above a certain temperature [6, 
100]. In order to investigate this abnormality, theoretical models and simulations have 
been developed that consider grain rotation [101], triple junctions drag [102, 103] or 
vacancy generation [104, 105] during the growth process. These models have been 
applied successfully to explain particular experimental observations; however, a full 
understanding of linear grain growth on the nanoscale still needs for further studies. 
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4.3 About vacancies 

4.3.1 Thermodynamic of vacancies 

The minimization of total energy implies that vacancies exist at thermodynamic equi¬ 
librium in any atomistic configuration. For the simplest case of pure metals, omitting 
the vacancy-vacancy interactions, the enthalpy and entropy variation due to adding 
vacancies will be, 

AH = X V AH V (4.1) 

and 


AS = X V AS V + R a (X v lnX v + (1 - A„)ln(l - X v )), (4.2) 

where X v is the molar fraction of vacancies, A H v is the increase in enthalpy and A S v 
is thermal increase in entropy per unit of mole of vacancies. At the equilibrium, 


dG 

dX, 


(4.3) 


where G = Gq + AH — TAS is the Gibbs energy function. Thus, one obtains the 
equilibrium concentration of vacancies as, 


X^ q = exp 



exp 


f -A H v \ 

V RgT J ' 


(4.4) 


The first term in right-hand side is a temperature-independent constant almost 
equal to 3 for metals. The second term, however, significantly increases with tem¬ 
perature and close to the melting temperature, vacancy concentration may rise up to 
10 -3 . In reality, however, the concentration of vacancies is usually non-equilibrium. 
This is due to the existence of defect sources and sinks which can produce vacancies 
or annihilate them. For instance, vacancies may be generated during removal of grain 
boundaries which will be discussed later in this chapter. Another possible source for 
extra vacancies are micro-voids which are formed during non-equilibrium production 
processes such as severe deformation or sintering. 


4.3.2 Strains related to the vacancies 

A vacancy maintains usually an empty volume which is smaller than that of a single 
atom and causes a small volumetric distortion in the lattice. Assuming vacancies to be 
the secondary element in the crystalline’s lattice and apply Vegard’s law, the elastic 
strain corresponding to the vacancies is given, 

el j = (4.5) 

ai 

with a v and a; being vacancy and substance lattice parameters respectively, c v the local 
vacancy concentration and the Kronecker symbol. In the presence of stress gradient, 
vacancies are able to move. This produces a micro-flow of atoms like it happens during 
creep. In general, vacancies influence the bond strength of materials and can change 
the elastic constants of materials as well [106]. In the following studies I focus on 
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the interdependency of vacancies and mechanical properties only through respective 
eigenstrain variations. 


4.3.3 Vacancy generation 

Vacancies are thermally enhanced defects. Their existence can be prescribed as natural 
response of lattice against increasing volume, decreasing pressure or increasing temper¬ 
ature. Vacancies may also be created due to phase transition or any other structural 
changes which are accompanied with local volume changes. Estrin et al. [107] have 
pointed that the reduction of grain boundary area during grain growth may cause con¬ 
siderable release of vacancies. They calculate the rate of total free energy change during 
this event, 


G 


3 a NkT 

V + -^{c-c eq )c, 


2 R 2 


c eq 


(4.6) 


where a is the interfacial free energy, R is the average grain size, V the interfacial 
velocity and c eq the equilibrium vacancy concentration. The first term is related to 
interfacial changes and the second one is the excess free energy rate due to excess 
vacancy concentration (c — c eq ). Grain growth continues only if G < 0. They have 
obtained the rate of concentration change as, 


C=^ V ~^-( c ~ c sink) (4-7) 

where e is relative excess volume of interfaces, i] is interfacial thickness, D v is the 
vacancy diffusivity, d is sink spacing and c s i n k is vacancy concentration at the sinks. 
They propose some other grain boundaries as th sinks. The first term in the right hand 
side express vacancy generation and the second term gives the vacancy removal rate. 
The solution for Equation 4.7 is approximated as, 


C = Cj + 



(4.8) 


where c* is initial vacancy concentration in bulk. Substituting this equation into Equa¬ 
tion 4.6 complete the calculations: 
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a 
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V + 24 



NkT 

c eq 


t < 0. 


(4.9) 


Inequality 4.9 says that vacancy generation increases the total free energy of the 
system. Therefore, above a certain time t* grain growth is stagnated. Thereafter, the 
system needs some time for redistributing the excess vacancies. In reality both of these 
phenomena happen simultaneously which results in new kinetics which is observed in 
related structures, e.g. Ns materials. 
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4.4 The model for interaction between interfaces and va¬ 
cancies 

I already mentioned that vacancies maintain a small amount of free volume and dis¬ 
tortion. Thus, the local density of atoms varies slightly from place to place within 
the substance. Assigning an unite volume to each atom it is possible to express the 
free volumes in terms of concentration. Considering a pure substance, one can define 
vacancy c v (x,t ) and atoms c a (x,t ) concentration fields with, 


c v (x,t) + C a (x,t) = c v0 (x,t) + c a o(x,t.) = 1. (4-10) 

c v o(x,t.) and c a o{x,t) are corresponding reference values 1 . Equation 4.5 maintains 
for the concentration field, 


e ij = o^_o, ^ 

ai 

On the other hand, interfaces are attributed by geometrically-required deformations 
(strain), corresponding to the mismatch between the grains. The amount of the strain 
is determined by local arrangement of the atoms, but it can be taken as a constant 
value in average. This deformation represents the free volume at the interface which 
can be significantly large for nanocrystalline materials as mentioned in Section 4.2. In 
the phase-field notation interfaces can be specified with a proper function h e ((f>) of the 
phase-field variables. Therefore, one can express distribution of geometrically-required 
strain as, 


ef = h e (</>) 4*S ij - (4-12) 

h e (<f>) is a non-zero interfacial function and e^* is a constant. The integration over 
Equation 4.12 will be equal to the e g (the total interfacial strain) which is set propor¬ 
tional to the total free volume of the interface. h e ((f>) = </>( 1 — (f)) is chosen. Thus, across 
a dual interface, 

e 9 = e t [ ‘K 1 - 4>) dx = ^e°*• (4.13) 

•7 — 00 * 

For multiphase contributions I use a linear summation of all dual contributions as, 

4 “ = < 4 - 14 > 

i¥=j 

Summing up above contributions, the total eigenstrain will be e'v + and the 
elastic energy density obtains 

fmech = I( e v _ ( e *J + e * g ij ))C ijkl ( C w - (e“ + ef 1 )). (4.15) 

T At constant volume, the density field p(x,t) = ^r- with a reference density po(x,t) (for thermo¬ 
dynamically stable bulk phase at a certain temperature) is an identical definition. Thereupon, one 
obtains p(x, t ) = c °(f = 1 ~ C ^V’^ where V a is the volume of a single atom. 
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This mechanical energy acts against interfacial motion and hinders the growth (See 
Section 4.5.2). At the equilibrium, however, the interfacial deformation can be com¬ 
pensated with the diffusion of vacancies which means + e*g l:) —> 0. If the mechanical 
distortion is fully compensated, the motion of the interface will be free from mechanical 
effects and only curvature-driven. 

Once the interface moves, however, the balance between the two sources of the 
deformation breaks: the geometrically-required deformation tends to move with the 
interface velocity while vacancies migrates via diffusion. 

From the equation of diffusion it is known that atoms always flow in a manner which 
would decrease the chemical potential gradient, 

j = M a V n, (4.16) 

where M a is the atomistic mobility and /i = ^ is the Gibbs chemical potential. This 
equation implies that the flux goes to zero if the chemical potential becomes equal 
everywhere. Here F is the total free energy functional which the elastic energy is a 
subset of that. Substituting Equation 4.15 into Equation 4.16 one obtains the flux of 
vacancies driven by elastic energy heterogeneities, 


j = M a V = M a V (^L_^i § ij C ijkl (e kl - (e kl + e* kl ))\ . (4.17) 

This equation demonstrates that the gradient of hydrostatic stress produces a 
mechanochemical potential which drives vacancies. In other words, vacancies will mi¬ 
grate in the stressed region around a moving interface seeking mechanical relaxation 

around the interface — —* 0). In the simulations, I treat density variation with 
a concentration field c v (x,t), which express the local free volume, and solve its time- 
variation, 


c = V.j = M a 


a r ai 


V.V 


ai 




(4.18) 


The constant M a is a temperature-dependent quantity equal to atomistic mobility 
in the bulk material. The Vegard’s factor av ~ ai may also depend on the temperature, 
but it is taken as a constant in the current model. 


4.4.1 Mechanical solution for interfacial deformations 

Existence of the strain in the interface causes elastic stress in the body of the system 
which can influence different events, such as mass transport phenomena, in a long- 
range distance within the bulk of the grains. Hence it is important to know about the 
distribution and amplitude of the stress within the local volume around the interface. 
For the sake of feasibility for an analytical approach, a single spherical grain embedded 
in an infinite substrate is considered which is separated from the substrate through 
a deformed interface. The analytical solutions, given in Appendix B, are plotted in 
Figure 4.1. Here it is assumed that materials inside, within, and outside the interface 
are isotropic with the same elastic constants. The interface is expanded volumetrically 
by 1 % ( e* B rr = e*£ 6 = ef* = 0.01). 
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Figure 4.1: Analytical solution for interfacial distortion are shown. Volume B (0.5 /jm 
width and 1 % volumetric expansion) is the interface between grain A (1.5 /iin radius) 
and substrate C. All phases are assumed to have same isotropic elastic constants (Young 
modulus of 208 GPa). The stress distribution is shown in the graph (right). It is found 
that there is no stress build in the grain interior due to the interfacial distortion. 


Surprisingly, it is found that regardless of its elastic constants, the grain does not 
feel any stress in this condition. The interface volume is under pressure and the solution 
within the substrate is similar to the Eshelby solutions (see Section 4.5.1). In other 
words, the volumetric deformation of the interface does not cause pressure within the 
grain with isotropic elastic properties. 

If the elastic constants of the interface and the substrate are different from each 
other ( Eb ^ Ec), the grain (domain A) will also be under stress (Figure 4.2). The 
source of the pressure inside the grain, indeed, is not the interfacial deformation, but the 
heterogeneous elastic properties of its surrounding medium. In general, if any obstacle 
in the interface or substrate interfere to the stress fields, the grain’s interior will feel 
the pressure as well. In the polycrystalline structures, each grain is surrounded with 
many other grains and grain boundaries which provide the heterogeneous media and 
result internal stresses inside the grains. It is worth noting that the sign of the stress 
depends on the sign of the eigenstrains, whether there is expansion of contraction at the 
interface. The motion of the interfaces is also a source of heterogeneity that significantly 
contributes to the stress distribution which will be discussed in Section 4.6.1. 


4.5 Mechanical contribution in the diffuse interface evo¬ 
lution 

In the MPF model, the volume of simulation is considered to be at mechanical equilib¬ 
rium which assumes mechanical stresses propagate instantaneously through the system. 
In the following, brief details about the mechanical equilibrium, the mechanical solver, 
its verifications, and the mechanical driving force related to the model are introduced. 
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Figure 4.2: Analytical solutions within and around a spherical grain (1 % volumetric 
expansion at the interface), for softer substrate, Eb > Eg, (left), and softer interface, 
Eg > Eb, (right), are shown. The elastic property of the spherical grain do not influence 
the results. 


4.5.1 Eshelby test of spectral mechanical solver 

Since an Eulerian method for finite difference volumes is used, it naturally does not 
allow for expansion or contraction of the simulation box. To get ride of this artificial 
effect, the volumetric changes of the box should be taken into account. This is done 
by exploiting an iterative Fast-Fourier Transformation (FFT) method which mimics 
the volume changes of the simulation box: first the strains are calculated for fixed 
system, thereafter one subtracts/adds an average volumetric strains (diagonal elements) 
from/to all discrete volumes iteratively until the total hydrostatic stress over whole box 
vanishes. This reflects a free boundary system allowed to deform at its surface or equally 
an infinite medium without stress/strain at far distance. Without this scheme the finite 
size effect will strongly influence the mechanical equilibrium. 

In order to verify applied homogenization scheme over the diffuse interface, as well 
as the free volume expansion scheme mentioned above, the analytical Eshelby solution 
for a spherical inclusion of isotropic elasticity to which it is possible to compare the 
simulation results is considered. Figure 4.3 presents the cr 11 and a 33 in a plane normal 
to 1 axis (x) for an inclusion of radius 17 nm with 1 % volumetric expansion (e* = 0.01 
along each axis). Both the inclusion and matrix are isotropic with C 11 = 289.8 GPa, 
C 12 = 163.4 GPa and v = 0.33 (as reported for pure cobalt). The simulation box is 
129 nm on each side and the interface width is 5 nm. Analytical solutions tell us there 
is a hydrostatic pressure inside the inclusion with constant diagonal elements of, 

a™ = -h*{C 11 + 2C l2 )^f^~ r<R p . (4.19) 

O -L IS 

R p is the inclusion radius. Within the matrix and in the chosen plane (perpendicular 
to axis 1) the stress components are, 


er 


li 


= -2cj 22 = —2<r 33 


-\e* {C 11 + 2C 12 ) 
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1 - 2v 
1 - v 



r > R p . (4.20) 
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Figure 4.3: Stress distribution around a spherical inclusion with 1 % volumetric expan¬ 
sion. 


Outside the interface region and inside the grain, the numerical solutions converge 
pretty well to the corresponding analytic solutions. Deviation at the interface strongly 
depends on interfacial width with respect to the size of the inclusion. The tangential 
stress components are less accurate. The deviation for diagonal elements, however, is 
still acceptable for our phase-field calculations. 

4.5.2 Mechanical driving force 

In the phase-field modelling, the contribution of non-interfacial energy terms are sep¬ 
arately calculated. This is already discussed in Chapter 1. The general derivation 
applies for mechanical contributions, 

^=- Uk - ^) ,mech ’ (4 - 2i) 

where depending on homogenization scheme, the driving force would be slightly dif¬ 
ferent (See Appendix C). The new model, however is slightly different such that the 
interfacial area has its own mechanical property. However, this does not influence 
the mechanical nature of the system. Assuming continuity of stress and displacement 
all over the space and homogeneity of elasticity, the Voigt homogenization scheme is 
applied. Therefore, 

AG™p ch = - e ° g *6 ij C ijkl (e kl - (e kl + ef))^ - <j> p ). (4.22) 

For the current case, the geometrically-required interfacial eigenstrain is distributed 
through the diffuse interface such that it has a maximum at the centre and vanishes 
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inside the phases which is adopted to the physical description of a symmetric interface. 
As a result, the mechanical driving force across the interface shapes with two extrema. 
The absolute values for the driving force overall is relatively small comparing to the 
extrema. This produces instabilities at the interface. In order to avoid this, a proper 
averaging procedure along normal direction of the interface is developed. The details 
of the averaging method is described in Appendix D. 


4.6 Demonstration of the model in the MPF scheme 

4.6.1 Shrinkage of a single nanograin 

In order to demonstrate the new model introduced in Section 4.4, a single spherical 
grain embedded into a substrate has been again considered. Having the same mate¬ 
rial constants everywhere, the interface maintains geometrically-required eigenstrain 
related to the free volumes and the conventional interfacial energy. The simulation 
starts from equilibrium, where the geometrically-required eigenstrain at the interface is 
compensated with the vacancies. The isotropic elastic constants of pure cobalt are 
taken as C 11 = C 22 = C 33 = 289.8 GPa, C 12 = C 13 = C 23 = 163.4 GPa and 
C 44 = C* 55 = C 66 = 63.2 GPa. The initial concentration at the bulk is 2 x 10 -4 . 
The interface energy is chosen to be 1 Jm“ 2 and the interfacial width is 6Ax where 
Ax = 0.2 nm. The mobility of vacancies is 10 _25 m _2 J _1 . Each vacancy is assumed to 
have 0.9 of the size of an atom, then the Vegard’s coefficient will be -0.1. The interfacial 
mobility, on the other hand, is fixed to ICG 19 m 4 J _1 s _1 . The numerical time step is 
0.1 s. 

The spectral mechanical solver obtains the mechanical equilibrium and the box ex¬ 
pansion as described in Section 4.5.1, in each time step. The condition for relaxation 
is if the absolute value of hydrostatic pressure is less than 10 3 Pa or the hydrostatic 
strain is less than 10“ 6 . Then, the diagonal components of stress matrix are embed¬ 
ded into Equation 4.18 and gave the hydrostatic stress gradients and the subsequent 
concentration increment in that time step. 

Figure 4.4 shows the results of this simulation. Initially, the mechanical distortion 
at the interface are assumed to be fully compensated by diffusin of atoms and vacancies. 
Here, it is needed to introduce the concept of “mobile atoms” which are in contrast to 
the vacancies. In other word, where the vacancies diffuse out from, the mobile atoms 
diffuse in. The number of mobile atoms are equal to the vacancies. But there are many 
potential mobile atoms, as each vacancy in surrounded by many atoms. We might use 
this concept in our future discussion. 

Indeed, depending on the sign of the geometrically-required eigenstrain at the in¬ 
terface, vacancies may diffuse into or out from the interface. It is important to be 
understood that the local distribution of stresses can be both positive or negative, re¬ 
gardless of positive eigenstrains at the interface. Since a positive build-in eigenstrain 
is assigned to the interface, vacancies are needed to diffuse into the interface to com¬ 
pensate them. The sign of the stress, strain and concentration components come out 
from derivations discussed in the model. The negative sign of the concentration field 
mentions the absent of vacancies and means a compressed structure (higher density), 
while the positive value presents the vacancies themselves. 

As it is expected the vacancies follow the stress gradient within and around the 
interface (Figure 4.4). Regarding the direction of the motion, geometrically-required 
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Figure 4.4: Evolution of concentration and pressure fields for a single spherical grain. In 
front of the moving interface there is a compressed area where the vacancies are adsorbed, 
while behind the interface, some vacancies are left. 
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eigenstrains will be annihilated behind the interface while should be build in the front. 
Thus, in the front of the moving interface vacancies are absorbed to the interface 
and behind the interface, mobile atoms migrate to compensate the stress. In other 
words, behind the interface vacancies are pushed to the substrate. Due to shortage 
of simulation time (small grain) and changing-curvature which varies the total driving 
force, it is not possible to observe a distinguishable stable motion, where the stress and 
concentration fields are in balance. The pressure difference inside and outside the grain 
has a gradient and increases for a shrinking grain. This can be understood also in the 
frame of local non-equilibrium density difference between the two sides of the interface. 

The kinetics of shrinkage is shown in figure 4.5 for different mobility factors. If 
there would be no diffusion at all, it is seen that the growth is slowed down with the 
mechanical deformation at the interface which is linear in the period of simulation. In 
contrast, when the diffusion is considered, the shrinking accelerates. If the mobility of 
atoms increases, which is possible by means of temperature, the compensation happens 
faster and eases the motion of the interface. For high enough mobilities the growth 
happens just same as a normal shrinkage due to the curvature with a parabolic kinetics. 

-9 



Figure 4.5: The kinetics of shrinkage of a single spherical grain is compared for three dif¬ 
ferent atomic mobilities with Ml= 10 -25 m -2 J -1 s -1 with no diffusion and normal shrink¬ 
age cases. At high mobilities the kinetics of shrinkage accelerates and get close to normal 
shrinkage process. 
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4.7 Simulation of polycrystalline nanostructured materi¬ 
als 

The general procedures of the current simulations are similar to what were previously 
described in Section 2.5.1. Due to computation costs for mechanical and diffusion 
solvers the simulation box is limited to 200 grids in each side with a grid spacing of 
0.2 nm, and interfacial width of 1.2 nm. A periodic boundary condition is used for the 
simulation box containing 1000 initial grains with almost the same size. The interfaces 
were assumed initially to be at mechanical relaxation. The mobility of vacancies M a is 
subject to change which reflects the annealing temperature, but the interfacial mobility 
is kept constant. Other parameters are same as it is described in Section 4.6.1. 

4.7.1 Simulation results and discussions 

4.7.1.1 Kinetics of nanogrowth 

The interfacial energy and the mechanical energy due to the deformation at the inter¬ 
faces directly contribute into dynamics of the interface. But the diffusion process which 
influences the total eigenstrains within both the interface and the bulk contributes into 
the kinetics of the grain growth, indirectly. Figure 4.6 presents the evaluation of av¬ 
erage grain size in different systems. The results are presented for the interfaces with 
30 % of free volumes. Turning off the diffusion (very low temperature), it can be seen 
that the growth almost stops. In this case, once the interface wants to move (together 
with the geometrically-required eigenstrain) a large amount of mechanical energy pills 
up which hinders the growth. This is where the low speed growth or no growth at all 
is expected. In this case the interface motion is partially or fully locked by mechanical 
distortions. 

For the atomic mobility a reference value M a o = 10 -25 m -2 J -1 s -1 is assumed. 
Simulation results for M a o, 10M a o and 100M a o are given in Figure 4.6. Once the 
diffusion turns on, vacancies move at the stress gradients, the eigenstrains compensate 
and the mechanical distortion at the interfaces can relax. This let the interfaces to 
move and therefore the growth happens. At this state, the motion of the interface is 
controlled through both the curvature of the interface and the diffusion of the vacancies. 
The speed of the interfaces depends on the rate of the diffusion: if the diffusion is 
slow the bottleneck of the growth kinetics will be vacancy diffusion and the growth 
is diffusion-controlled. When the atomic mobility is equal to 100M a o the mechanical 
distortions fully compensate at the interfaces and the kinetics of growth almost match 
with the normal grain growth, which means the growth is again curvature-driven. 

The importance of free volume redistribution during grain growth is already men¬ 
tioned in Sections 4.2 and 4.3. Estrin, Gottstein and Shvindlerman [104, 105, 107] 
suggested that during the grain growth interfaces release vacancies into the bulk that 
increases the total free energy of the system and can prevent the growth. Generated va¬ 
cancies were supposed to diffuse through the bulk to the possible sinks which should be 
another interfaces. This model has been applied to the grain growth in pure nanocrys¬ 
talline [7] and thin films [108]. The present work, however, considers another aspect of 
vacancies which are influenced by the elastic energy potential accumulated at the in¬ 
terface. Shvindlerman [109] and Faulkner explain the thermodynamic of the vacancies 
at the interfaces where the chemical potential of interface can reduce by adsorption 
of the vacancies. This is very similar to the concept presented here, except that the 
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4.7 Simulation of polycrystalline nanostructured materials 



Figure 4.6: Nanocrystallined cobalt are simulated and compared to the normal grain 
growth of the same system. The grains stop growing if there is no diffusion. For different 
diffusion rates the kinetics of the growth changes from the diffusion-controlled process to 
a curvature-driven motion of the interface. The interfaces are considered to have 30 % of 
free volumes. 
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4. VACANCY DRAG: A MODEL FOR KINETICS OF THE 
INTERFACE ON THE NANOSCALE 
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Figure 4.7: The average grain size of pure cobalt nanocrystalline at the end of 1 hour 
annealing process. The growth mechanism changes at a critical temperature 803 K. Ex¬ 
perimental results by Song et al. [6]. 


mechanical subset of the chemical potential at the interfaces is considered to reduce 
by gathering of vacancies at deformed interfaces. One advantage of this model over its 
thermodynamic alternative is that the long-range interaction between interfaces and 
mobile defects is explained. 

In the annealing experiments of nanocrystallined pure cobalt reported by Song et ah, 
the transition from the linear to normal grain growth is observed starting at 803 K is. 
Estrin et al. [107] suggested a temperature-dependent critical grain size of about 200 nm 
above which the growth is no more inhibited by vacancy generation and normal grain 
growth occurs. In the MPF model, although the direct simulation of the transition from 
linear to normal grain growth is difficult due to abundant computational requirements, 
the nature of the transition can be discussed as follows. For very small grains, the 
regular curvature-driven motion would be faster than the maximum velocity allowed 
by vacancy diffusion. Hence, in this regime, the diffusion-limited motion of interface 
causes the linear grain growth kinetics. Indeed, the vacancies which move slowly at 
low temperature impose an indirect drag force on the interface and slow it down. 
Once the grains reach a critical size above which the curvature contribution to the 
motion of the interface is lower than the diffusion-controlled kinetics, the controlling 
mechanism changes. Obviously, the transition size is highly temperature-dependent. 
For a constant annealing time the grains may not reach the critical grain size at low 
temperatures which means the growth process is only a linear growth mechanism. In 
contrast, at high temperatures a fully normal grain growth is expected. Between these 
two extreme cases, there is a transition range where both linear and normal regimes 
contribute to the growth process. 

The discontinuity on the nanoscale grain growth observed experimentally for pure 
Co between 800 to 900 K. Figure 4.7 presents the exprimental results obtained by Song 
et al [6]. The activation energy for migration of vacancies is similar to the self-diffusion 
activation energy in the bulk of the materials. The self-diffusion in cobalt is expressed 
in following Arrhenius equation with a large range of temperature [110], 
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D Co = 1-287 x 10 _4 exp 


-296542.9 

R 9 t 


m 2 s 4 . 


(4.23) 


Applying Einstein relation -which neglects double-vacancy diffusion and assumes a 
dilute concentration of vacancies- one obtains, 


Ma = M Co = (4.24) 

Thus the corresponding temperatures for M a o, 10M a o and 100M a o are « 898, ~ 957 
and ~ 1095 K which are above experimental values. Although the grain size in our 
simulations are much smaller than reported values in Figure 4.7, it still provides a 
fair explanation of the possible mechanisms of transition from low-speed grain growth 
into normal grain growth. From Molecular Dynamics (MD) calculations for copper 
which were performed parallel to this study in our group [111] we observed that the 
activation energy of the shrinkage of a single spherical nanograin is about 10% of the 
activation energy for self-diffusion, this is also reported in some experimental researches 
[112, 113, 114]. This low activation energy suggests the interfacial diffusion rather 
than bulk diffusion. In this situation, the speed of the interface is determined by the 
frequency of atomic jumps across the interfaces. Since we do not assign any special 
kinetic aspect to the atoms at the interfacial region, the activation energy is everywhere 
the same which can be one reason for deviation form experimental observations. The 
MD simulations, however, confirm that the moving interface may leave some vacancies 
and stacking-fault area (in the case of Cu) behind, which affirm the redistribution of 
the free volumes. 

4.7.1.2 Microstructure and grain size distribution 

Figure 4.8 compares the results of two simulation boxes with 10M a o and 100M a o atomic 
mobilities. Both simulations start from the same initial structure. Later on, not only 
they obtain diverse growth kinetics (as shown in figure 4.6), but also the trajectory of 
microstructural evolution are different. This can be understood because of long-range 
interactions between diffusion and stress fields. In fact, the flux of atoms in each point 
of the space may be influenced by more than one interface. 

Figure 4.9 shows the concentration and pressure fields at a cross section view for 
the simulation box with 10M a o. The location of the interfaces can be easily recognized 
where the concentration field reaches its minima. Clearly, the concentration of vacancy 
at the interface is very high. At the triple junctions the concentration is even higher, 
proportional to the amplitude of deformation. The negative sign for the concentration 
shows the absent of the vacancies or presence of mobile atoms. This can be interpreted 
as densification due to pressure. At the beginning, the pressure distribution is very 
heterogeneous and mostly localized at the interfaces and its adjacent regions. Later on, 
the interior part of grains also feels more of a pressure, while the interface volumes are 
highly compressed. The importance conclusion is that the pressure is not a direct result 
of interfacail distortion, but a consequence of migration of the deformed interfaces. 

The experimental measurements show that the interfacial stresses at the nanoscale 
structure can easily exceed lGPa [115]. At the triple junctions one can see that this 
value can be remarkably larger (up to 40 times). Of course this occurs in a very short 
time, when the interfaces suddenly move and the instant mechanical deformation is 
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4. VACANCY DRAG: A MODEL FOR KINETICS OF THE 
INTERFACE ON THE NANOSCALE 



Figure 4.8: Cross section of simulation boxes for 10M a o (left column) and 100M a o (right 
column) atomic mobilities. The images in each row belong to the same time step (from 
top to below initial, 100s, 200s and 300s). The grains “1” and “2” can be clearly compared 
in two set of simulations. 
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Figure 4.9: Concentration and pressure fields are compared for 10M a o after 10 (left) and 
1000s (right). 
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<R> ~ 2.1 nm <R> ~ 3.1 nm <R> ~ 3.9 nm 



Figure 4.10: Size distribution for normal grain growth (top row) and diffusion-controlled 
grain growth (below row) at three different steps. For comparison, the distributions are 
presented for almost the same average grain sizes. 


build up. Thereafter the amplitude of pressure decreases. Vacancies, however, can not 
catch the large stress gradients unless diffusion happens fast. At the later stages of the 
growth, when the diffusion and stress are almost balanced and their changes are small, 
the pressure is reduced and nearly equal over the system. In this step, the concentration 
is localized at the interfaces to relax the distortions as much as possible. 

The grain size distributions were compared for diffusion-controlled regime in Figure 
4.10. For better comparison, the distributions are presented for almost the same average 
grain sizes. Since statistics are pretty poor (limited number of grains), it is difficult to 
investigate whether transition to a self-similar regime for nanograin growth happens or 
how it may differ from normal regime of growth. However, the evolution of diffusion- 
controlled growth is distinguishable in general. At the nano-regime, tendency towards 
a rather narrow distribution is observed: population of the (relatively) small grains 
rapidly decline and a plateau forms in the left side of the graphs while the relative 
number of grains at the average-size remains intense. It seems that below a certain 
size, the curvature overcomes the mechanical barrier for motion. Another possibility is 
that the stress fields may neutralize each other through long-range interaction which is 
doable across the small grains. 


4.8 Chapter conclusion 

A phase-field model has been verified and applied for grain growth on the nanoscale 
in which the existence of extra interface volume possessing mechanical distortion is 
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4.8 Chapter conclusion 


taken into account. The interaction between interfacial stress and vacancies, as mo¬ 
bile defects, introduce a new regime for dynamics and kinetics of the interfaces. The 
geometrically-required deformation at the interface builds mechanochemical potential 
gradient, the compensation of which seeks vacancies to diffuse-in from the adjacent 
volumes. Therefore, the kinetics of interfacial migration found to be controlled by bal¬ 
ancing between diffusion of the vacancies, as well as the curvature of the interface. At 
low temperature, slow diffusion is the bottleneck of the grain growth which conducts 
low-speed linear growth regime on the nanoscale. At this stage, even the curvature 
driving force is substantial (small grain size) the motion of the interface is locked by 
mechanical distortion, the relaxation of which requires rapid diffusion. At high tem¬ 
peratures, however, the diffusion speeds up and meets the mechanical relaxation at 
interfaces. The growth, therefore, is curvature-driven in this condition. Between the 
two extreme cases, the balance between two controlling mechanisms determine the ki¬ 
netics of the growth. Once the mobile defects can not catch the moving interface, 
an internal pressure is build up within the grains and grain boundaries proportional 
to the local density variation. This pressure is long-range and can contribute to the 
neighbouring grains and interfaces. It is observed that in this situation relatively small 
nanograins are very unstable and grains with the average grain size persist more than 
any other grains. 
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Appendix A 

A Brief note on code 
developments 


To accomplish present thesis, much effort has been made to create a new efficient and 
quick program which enable us to simulate the physical problems of interest. This 
code has been developed in C++ environment and has a modular structure, in which 
distinct classes for phase-field, concentration and elasticity are defined. Each class is 
meant to have its own local variables/functions as well as access to needed global vari¬ 
able and functions, and to be stored in a separate hie. The classes also are coupled 
within specific functions to express special effects such as vacancy-stress relationship. 
The main() function then calls the initial functions (needed for reading inputs, ini¬ 
tialization, memory reservation,...) and maintains main-loop which handles simulation 
steps and writing outputs (see flowchart Figure A.l). The results of current code de¬ 
velopments are collected in OpenPhase [116], the open-source package for phase-field 
simulations. In the following I will give a brief explanation about the code structure 
and its developments. 


A.l Computation of derivatives 

In MPF equations, time and space derivatives of phase-fields , concentrations and dis¬ 
placements are required. These contributions may interact with each other, depending 
on the problem during the computations. For example, the normal vectors to the inter¬ 
faces computed in the phase-held space may be employed for any other events such as 
averaging of the driving forces, evaluation of the interface stiffness etc. Here I present 
few more details about derivatives. 

Phase-field: The gradients of phase-fields are calculated using their first nearest 
neighbours and for their Laplacian, I apply an isotropic stencil with include 26 neigh¬ 
bours (first nearest neighbours plus diagonals) [66]. For time derivatives, an Eulerian 
finite difference scheme for both phase- and concentration fields is applied. The phase- 
fields are summed up with all other driving forces created by mechanical and diffusion 
contributions and updated in each time step. Phase-field computation is limited to the 
interfaces and their first neighbours. 

Concentration field: For the concentration, a general diffusion function should be 
solved which requires chemical potentials as an input. In contrast to phase-fields, the 
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A.l Computation of derivatives 





Type definition 


main() function 





Phase-fiield increaments 



Figure A.l: Simple flowchart of computation algorithm. 
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A. A BRIEF NOTE ON CODE DEVELOPMENTS 


concentration should be solved over whole simulation box. The time step of diffusion 
steps is usually small and the bottleneck against speed of the calculations. 

Mechanical fields: Strains and stresses are the output of the mechanical solver. 
We use a spectral solver based on FFT method which is rather fast but limited to 
periodic boundary conditions. Its iterative scheme allows for tuning precision of the 
solution depending on the problem and the performance concerns. The grid points 
for mechanical solver match the phase-field grid points. Due to the symmetry, there 
are only 6 independent components for strain and stress at each point in the space, 
which is still a huge amount of data, derivatives of elastic energy can be input for both 
concentration and phase-fields classes. 


A.2 Memory structure 

One of the challenges for simulation studies presented here was the large amount of 
data which should be stored or written as output. This requires large amount of storage 
capacity as well as extra computation time for writing them. In the MPF model N 
phase-field variables are assigned to each point in the space requiring N x nX x nY x nZ 
units of the memory which is a large number in our simulations. Technically, however, 
the number of existing phases 1 at each point in space is much less than N and frequently 
there is only one phase (inside each domain). Therefore, we designed a dynamic storage 
algorithm which only stores the existing phase-fields but not all of them. This is 
handled using a C++ vector container which allows for easy resizing, looping, adding 
and deleting of each individual phase value(s), due to its pointing mechanism. A 
phase-field list contains C++ structures, each block having a phase number and its 
phase-field value. Thus, the size of the list is equal to the number of available phases. 
The schematic description is given in Figure A.2. In this way we succeed to increase 
the size of the simulation domain which was crucial for the studies of normal grain 
growth. 


A.3 Inputs and outputs 

All input parameters were read from an ASCII file which eases the interaction with 
the code. Each module has its own input file where only its relevant parameters are 
provided. The type of variables are mostly defined in a Settings class which take care 
of general input parameters. The phase-field variable is taken as float which is smaller 
than double (represented 32 bits against 64 bits) but still satisfies the expected accuracy. 
The concentrations, strains and stresses, however, should have a double precision. 

For initialization input of big simulation boxes, Voro++ program [67] was partially 
modified and used. The domains have usually sharp interfaces at the beginning which 
should be carefully treated at the initial step of simulation. For other cases, initializa¬ 
tions are written in individual predefined functions with diffuse interfaces. 

Writing output is one of the time-consuming steps of the simulations. Read¬ 
ing/writing functions have been developed such that they gather minimum amount 
of data. Phase-field values, concentrations, stresses and strains were the general out¬ 
puts of the simulations. In order to preserve accuracy, they were partially written/read 

1 In the phase-field simulations grains of the same material are distinguished by their local lattice 
orientation. 
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A.3 Inputs and outputs 



Figure A.2: Schematic demonstration of the storage structure of phase-fields. The phase- 
field value of phase p in the point (*, j. k) in the space is stored in a 3D vector of a structure. 
If the phase does not exist, the algorithm automatically returns the value 0. 


in binary format and stored in separate files and named according to their respective 
time-step. In this way, it is easy to continue interrupted simulations. Phase-field out¬ 
puts contain the number of phases and their values in all points of the domain. Outputs 
with “vtk” format were generated for visualization using Paraview software and the raw 
outputs feed the statistical analysis programs developed in MATLAB and C+-h 
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Appendix B 

Mechanical equilibrium around a 
deformed interface 


I obtain analytical solutions for mechanical equilibrium around a deformed interface. 
The simple model depicted in Figure 4.2 considers a single spherical grain of radius Ra 
separated with a diffuse interface of width Rb — Ra from an infinite substrate (domain 
C ). Deformation of the interface is assumed to be only volumetric and reflected in its 
eigenstrains. The properties of three domains are denoted with indexes A, B and C. 
Assuming continuity of stress and displacement, the analytical solutions in spherical 
coordinate (r, 9, 4>) are given by 
For r < Ra : 


a = 


E a 


(1 + u)(l — 2u) 
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For Ra < r < Rb : 
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B. MECHANICAL EQUILIBRIUM AROUND A DEFORMED 
INTERFACE 


where 
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The Young modulus for isotropic material is E = — — — = - — l 1 +y i = 

2(1 + v)C M . 
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Appendix C 

Homogenization methods at 
diffuse interfaces 


Since phase-field is about multi-phases (grains) in contact, it should be combined with 
proper homogenization methods to satisfy the continuity conditions for stress and dis¬ 
placement. A few well-known methods which are relevant to our studies are briefly 
introduced here. 


Voigt model In this model, it is assumed that the strain is uniform among all do¬ 
mains. This is true when two phases of similar stiffness are brought into contact. 
Applying this assumption into phase-field, the total strain at the diffuse interface con¬ 
sidered to be equal ( Attention: in this appendix I avoid the indexes for stress, strain 
and elastic constants to make the equations readable. Of course for demanding or 
unclear statements, I will explain them properly.): 

e = e a = e/ 3 , (0.1) 

and effective stress is expressed by weighting stresses of different phases, 

a = 4> a a a + <j)p<jp = (j) a C a {e - e*) + <t>pCp(e - e*p). (C.2) 

Therefore the effective elastic constants and eigenstrains will be, 

C e ff = 4>aC a + <t>pCp , (C.3) 

< ff = C~ f ) {4> a C a e* a + <j)pCpe*p) . (C.4) 

The elastic energy stored in the system also can be described in the terms of effective 
parameters, 


r^ = -(e-etff)C eff (e~e* eff ). 


(C.5) 
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C. HOMOGENIZATION METHODS AT DIFFUSE INTERFACES 


Reuss model If the phases in contact are very different in elastic moduli, they may 
expand or contract very differently against forces. In such case homogenizing of stresses 
within diffuse interface is rather useful, 


(7 — (J a — (Jp. 

(C.6) 

The strain can be averaged between two phases, 


£ = 4 e a + 0/3 e /3 = 4 > aC a 1 (7 + 1 £T. 

(C.7) 

The effective elastic constants and eigenstrains would be 



(C.8) 

and 


e eff e a 

(C.9) 


Substituting these formulations into Equation C.5 the elastic energy density obtains. 


Khachaturyan model In 1983 [35] Khachaturyan introduced a mixture model which 
takes weighted eigenstrain and elastic constants, 

e*ff = </> a e* + <p/3 ejg, C e ff = (j> a C a + (j)pCp. (C.10) 

The first assumption is similar to Reuss approach whereas the second assumption 
is a result of Voigt model. The stress is given by Hook’s law as, 

a = C e ff(e — = (4> a C a + <t>pCp)(e — <44 — <f>pe*p), (C.ll) 

and the elastic energy obtains again by inserting these materials into Equation C.5. In 
contrary to two standard models above, the assumptions of Khachaturyan model does 
not allow for definition of phase elastic energy separately [117]. 


Steinbach model Steinbach et al. [3] proposed another approach, applicable for 
their MPF concept. They assumed uniform stress across diffuse interface, but they 
consider the elastic strain to be linearly weighted with phase-field variables, 

o = a a = ap (C.12) 

and 


e — e* = 4>a(,e a ~ 4) + <f>p(e a ~ 

The effective elastic constants here are the same as 
proach, and the elastic energy is again as Equation C.5. 


4 )- (C.13) 

Equation C.7 from Reuss ap- 
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Appendix D 

Averaging driving force over the 
interface 


Discrete volumes in our numerical scheme may influence the accuracy of the computa¬ 
tions. Even sometimes, it may produce artefacts and instabilities at the interface. For 
instance, if there is multi-extrema function acting on the interface area, the calcula¬ 
tions may fail to catch some critical points or to miscalculate the overall contribution 
of driving forces. 

In the model, presented in Chapter 4, the driving force due to geometrically-required 
eigenstrains builds a two extrema function with large range but small average value. 
This is shown schematically in Figure D.l. Since the local values are so divers, the 
interface may be very unstable. In order to avoid these effects, I average the driving 
forces at the interface and redistribute them by a rather smooth function. 


+ 



Figure D.l: Schematic sketch of driving force function over an ID interface (see Section 
4.5.2). The total average value is much smaller than extrema. 

Averaging should be done in the normal direction over the surface or line of the 
interface defined as For averaging, theoretically one should integrate the values 

along the normal vector and divide by its length. In a discrete space, however, any 
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D. AVERAGING DRIVING FORCE OVER THE INTERFACE 



Figure D.2: 2D Schematic expression of the averaging. The rectangular (cylinder in 3D) 
over the interface (curved area) is parallel to the normal vector plotted for a certain point 
at the centre of the interface. 


defined vector may hit randomly some grid points. This is not numerically enough 
accurate, especially for averaging over heterogeneously distributed values. To resolve 
this issue, a certain distance around normal vector are attributed to it such that it covers 
all its nearest neighbours. This is similar to averaging over values inside a cylinder for 
which its axis is coincide with the normal vector. Figure D.2 gives graphical explanation 
for this method. Thus, the average value should be redistributed using a proper function 
over the interface. Two typical non-zero function are ^d(l — 0) or 1 — </>), the 

integration of which over the profile of an interface with DO potential is equal to 1. 
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